Stochastic equation for a jumping process with long-time correlations 
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A jumping process, defined in terms of jump size distribution and waiting time distribution, is 
presented. The jumping rate depends on the process value. The process, which is Markovian and 
stationary, relaxes to an equilibrium and is characterized by the power-law autocorrelation function. 
Therefore, it can serve as a model of the 1/f noise as well as a model of the stochastic force in the 
generalized Langevin equation. This equation is solved for the noise correlations ~ 1/t; the resulting 
velocity distribution has sharply falling tails. The system preserves the memory about the initial 
condition for a very long time. 
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I. INTRODUCTION 

A jumping process can be defined in terms of two prob- 
ability distributions which determine a jump size and a 
waiting time between consecutive jumps. One usually 
assumes that both distributions are independent of each 
other. Such process is often regarded as a generalized 
form of the random walk and used to describe a diffusive 
transport. That approach, known as the continuous-time 
random walk theory is able to account for various 
forms of diffusion, both normal and anomalous, by a suit- 
able choice of the probability distributions defining the 
processJ3- Power-law dependences are especially inter- 
esting 0, 111 ■ A stochastic trajectory characterized by 
jump sizes so distributed, exhibits a pattern typical for 
the Levy flights and features systems which reveal the en- 
hanced diffusion 0. On the other hand, long tails of the 
waiting time distribution (long rests) evoke the opposite 
effect: they are responsible for the subdiffusion 0|. 
Processes which possess such tails are often treated in 
terms of the fractional diffusion equation 0, H, S ^ • 

For uniform distribution of jumps in time, i.e. if 
the waiting time probability density has the exponen- 
tial form, the jumping process relaxes to some stationary 
equilibrium. The kangaroo process provides a sim- 
ple and well-known example. Instead of the jump size 
distribution, this process assumes a probability distribu- 
tion of the process value after the jump and, in addi- 
tion, a jumping rate which depends on the process value. 
An advantage of the kangaroo process from the point of 
view of possible applications stems from the fact that 
it can be easily constructed for arbitrary correlations. 
A need for models of correlated noises is obvious. For 
example, the long correlations, both in space and time, 
arise as a result of the fast modes removal procedure 
0,0,0]. Long tails of the correlation function emerge 
also in the relaxation process of a system coupled to a 
fractal heat bath via a random matrix interaction |l5j| . 
In those cases the stochastic dynamics obeys the general- 
ized, non-Markovian, Langevin equation and the Monte 
Carlo simulation of solutions requires a specific model 
of the noise. Unfortunately, the kangaroo process is not 
suitable to model noises with power law correlations: the 



distribution of the stochastic variable during the trajec- 
tory evolution is biased because the waiting time dis- 
tribution changes its shape when it is inserted into the 
generalized Langevin equation As a result, the re- 
laxation to the thermal equilibrium cannot be achieved. 

In this paper we consider a simple power-law correlated 
jumping process which is exempt of that difficulty. Our 
objective is not only to analyze the master equation for 
that process but above all to obtain the stochastic vari- 
able itself by solving a stochastic equation. Therefore the 
presented procedure can be utilized as a noise model for 
numerical simulations of the stochastic trajectories in the 
framework of the Langevin formalism. We define the pro- 
cess and discuss corresponding equations in Sec.IL The 
expression for the autocorrelation function is derived in 
Sec. Ill, whereas Sec. IV is devoted to the application of 
the process as a model of some specific form of the cor- 
related noise in the generalized Langevin equation. The 
main results are summarized in Sec. IV. 



II. DEFINITION OF THE PROCESS 

We assume that a stochastic process is step-wise, i.e. 
a process value x is constant within the time intervals 
{U,ti+i): x(i) = Xi for t G {U,ti+i). Jumping times 
are randomly distributed and jumping rate i^(x) depends 
on the process value. The size of the jump, defined as 
the difference between the values of x after and before the 
jump, is determined from a given probability distribution 
(5(i5x). Then the stochastic trajectory x(t) obeys the 
following equation 



Xi+l = Xi 



(1) 



where the waiting time r = i^+i — is governed by the 
Poissonian distribution: 



Pp(r) = ;/(x)e- 



(2) 



which determines the probability density that a jump 
occurs in the interval (r, r + dr). The initial condition for 
the Eq.(^|, x(to) = xq, follows from a given probability 
distribution Po(x). The Eq.lJ^I is stochastic because it 
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determines the time evolution of the stochastic variable 
X, in contrast to the master equation which can give us 
only probability distributions. The trajectory x{t) can be 
constructed step by step by sampling consecutive values 
of T and (5x from the distributions Pp and Q, respectively. 

The process is Markovian and stationary. The transi- 
tion probability ptrdii. that the process value is between 
X and X + dx at an infinitesimal time Ai, providing it 
was equal x' at t = 0, is given by: 



Pt,(x, A^|,x',0) = {1 - i^(x')A0<5(x' - x)- 
+iy(x')Atg(x - x'). 



(3) 



In the above expression we have utilized the fact that 
Ptr may depend only on time differences. The first term 
on the right-hand side of Eq.Q is the probability that 
no jump occurred in the time interval (0, At). The term 
i'(x')At means the probability that one jump occurred. 
The master equation for a probability density p(x, t) 
can be obtained by calculating the time derivative from 
p(x, t) and taking into account all possible initial values 



d_ 

dt 



p(x, t) 



- lim 

At->0+ 



j Ptrix,At I x',0)p(x',t)dx' -p(x,t) 



/At. 



(4) 



We get the master equation in the following form 
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p(x, t) = -i^(x)p(x, t)+ I Q(x - xXx')p(x', i)dx'. 

(5) 

The jumping process described above is still too gen- 
eral and thus we introduce additional restrictions. Let 
X be a two-dimensional vector, x = (xi,X2), with the 
unit length: |x| = 1. Therefore we require the norm to 
be preserved during the jumps. With these assumptions, 
the process can be described in terms of a single angle 
variable (/>: xi = cos((/>) and X2 = sin(0). For the proba- 
bility density Q we take the Gaussian: 



g(<5x) 



-(x- 



Ne' 



(6) 



where cr is a given width and the normalization con- 
stant N = l/[27r/o(l/cr^)] contains the modified Bessel 
function. We will demonstrate in the following that the 
jumping rate v determines the shape of the process au- 
tocorrelation function. It appears algebraic if v is given 
by the expression 



|sin(0)|" 



1 — a I cos((/))| 



(7) 



where < a < 1. Taking into account the above as- 
sumptions, we obtain the master equation Eq.© in the 
one-dimensional form: 



^p((^,<) = -K0M</',t)^ 



(8) 

The equilibrium solution of Eq.JSJ, P{4>), has to satisfy 
the condition i'{(j))P{(j)) =const. Therefore, ^(0) becomes 
quite simple: 



Since for the jumping rate Q) Jq^ l/i 
properly normalized. 

Numerical simulation of stochastic trajectories requires 
random numbers distributed like Q(5x), according to 
Eq.@. For that purpose we apply the rejection method 
which allows us to avoid evaluating complicated integrals. 
The algorithm is the following. First we sample uni- 
formly distributed random numbers 5(j) = 4> — (f)' from 
the interval (0,27r). Then q = Q{S(f)) is calculated and 
this value is compared with another random number, rg, 
uniformly distributed within the interval {QmimQmax) 
where Qmin and Qmax denote the minimum and max- 
imum values of Q, respectively, li rg > q then Scf) is 
accepted, otherwise it is rejected and the sampling pro- 
cedure is repeated. 



(9) 



1, P(0) is 



III. AUTOCORRELATION FUNCTION FOR 
THE JUMPING PROCESS 



The autocorrelation function of the process (ACF), 
C{t) = (x(O)x(t)), where the average is taken over the 
stationary distribution P(x), can be evaluated from the 
following expression |0 



C(t) 



x'(to)x(to +t)P(x, t\x')p{x', to)dxdx'. (10) 



The conditional probability of passing from x' to x dur- 
ing the time t, P(x, i|x'), can be obtained by taking 
into account all possible paths leading from x' to x and 
summing over the jumps p^ . The final formula for the 
Laplace transform of ACF can be expressed by the fol- 
lowing series 
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Cis) 



27r 



1^(0) 1/(0) +s Jo 1^(0) +s 1/(00+5 



oo „27r 



E 



COS(0 - 0o) Q{4> - 0fc_] 

i^(0o) + s v{(j)) + s 



—j2 ^ , ^ i^{(Pi~i )d(pi-i 



.1=2 



J/(0i_l) + S 



(11) 



Inverting the Laplace transform we obtain the final expression for ACF: 



-v{4,)t 

m =4 / ^-^d0+ 



-8N 



7r/2 /.Tr/2 
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' gCos(</)-0')/o' _g-cos 



) cos(0 - 00 " W + 

y 1/(0) - z/(0O 



(12) 



We are interested in the asymptotic bahaviour of C{t) 
for large t. In this limit the first term of Ea. ()12|l can be 
estimated easily. Because of the exponential dependence 
of the integrand on t, only a vicinity of = contribute 
to the integral: i/ < 1/t. Therefore the first term can 
be approximated by the integral exp(— 0"t)/0 c?0 ^ 
^1-1/q ^]-^g second term we first calculate the integral 

over 0: Jq^^[Jq^^ d(j)]d4>' . If we take the limit of large t 
in the inner integral, the exponential containing 0' can 
be dropped. Moreover i/(0) becomes negligible, compar- 
ing to 1/(00, as well as in the arguments of the cosine 
function. Then for any 0' > we have: 



-1/(0')* 



r.(0) - 1/(00 
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and the integral over can be easily evaluated. The re- 
quired time dependence is of the form which means 
that the second term falls with time faster than the first 
one. The same conclusion refers to the higher terms. The 
second term has a simple asymptotic dependence also on 
the kernel width a. Expanding the exponential functions 
over 1/cr and taking into account that Hm^^o Io{x) = 1, 
we find that the second term decreases like 1 /cr^ for large 
a. We finally conclude that the ACF can be well approx- 
imated by the first term of the Ea. l|12|) and its tail is 
algebraic: 



C{t) ~ t^-^/°' for t 



(13) 



Fig.l presents ACF for a — 0.5; C{t) was calculated 
from the definition, by means of single trajectories evo- 
lution according to Eq.Q, for a = 1 and a = 2.5. The 
equilibrium probability distribution ^(0) was taken as 
the initial condition. The result for the larger value of a 
agrees very well with the first term in Eq. H12|) and it can 
be parameterized by the function 



-8t 



8t 



(14) 




FIG. 1: Autocorrelation function C{t) for the jumping rate 
v given by Eq.JTJ with a — 0.5. Numerical simulations have 
been performed for a — 2.5 (solid line) and a = 1 (dotted 
line). The first term in Eg. 11211 is also shown (dashed line), 
as well as the parameterization by Eg. 1141 (dashed-dotted 
line) . 



The existence of long tails of ACF means that the 
power spectrum of the process, defined by the Fourier 
transform J-'(x) as S{lu) — |J^(x)p, is strongly enhanced 
at Lo = 0. The power spectrum can be obtained directly 
from C(t) by using the Wiener-Khinchin theorem jl7l |: 
S{uj) = T{C{t)). For 0.5 < a < 1 we get the following 
result: 



-l/a 



(15) 



Then our jumping process is characterized by the al- 
gebraic power spectrum and becomes the 1// noise for 
a — > 1. Overpopulation of small frequency values is due 
to the fact that the process is dominated by long waiting 
times between consecutive jumps. Such long intervals 
correspond to small values of 0, i.e. to the evolution 
along the Xi axis. A quantity s = 1/j/, which means the 
average of the Poissonian distribution ||2Jl, is well suited 
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to characterize long rests. The statistics of s is directly 
connected with the process value probability distribution 
P{4>) and, in accordance with that, the density distribu- 
tion of s in the equilibrium, "(/"(s), can be derived from the 
equation \'ijj{s)ds\ — In the limit of large s we 

obtain -0(5) ^ s~^/" and this result means that the Pois- 
sonian waiting time distribution with variable jumping 
rate can possess, effectively, power law tails. 

The jumping process with a — 0.5 resembles a deter- 
ministic dynamical system: the Lorentz gas of period- 
ically distributed hard disks. In this lattice a particle 
is elastically reflected by the discs and wanders freely 
among them. Free paths of the particle are infinite at 
directions parallel to the symmetry axes. The system 
is characterized by the velocity autocorrelation function 
with the tail 1/t, analogously to Ea. H14|) . and by the 
power spectrum S{lu) ^ |ln(w)|. However, the long free 
path distribution falls faster then its stochastic counter- 
part, hke s^^ [19J. 



IV. APPLICATION TO THE GENERALIZED 
LANGEVIN EQUATION 

If a Brownian particle is driven by a stochastic force 
with a finite correlation time, the time evolution of the 
velocity obeys the generalized Langevin equation .20., ,21i |: 



m- 



dv(i) dV{r) 



dt 



dr 



m / K{t-T)v{T)dT + F{t) (16) 



where V{r) is a position-dependent external potential, 
F{t) is a stochastic force and m denotes the mass of 
the particle. The integro-differential equation H16() can 
be solved numerically for any V{r) if we apply a con- 
crete model for the noise F{t). In the case V{r) = the 
Ea. (|16|l is manageable by Laplace transforms. We obtain 
the following solution: 



V(t) = i?(i)v(0) 



R{t-T)F{T)dT, (17) 



where the Laplace transform of the resolvent R{t) is given 
by the equation 



(18) 



In the Eq. (|16|l the usual damping term - proportional 
to the velocity - which appears in the ordinary Langevin 
equation, has been substituted by the retarded friction in 
the form of a memory kernel to ensure proper character- 
istics of the equilibrium, namely the equation should sat- 
isfy the second fluctuation-dissipation theorem [25| . The 
kernel K{t) has to be proportional to the noise ACF C{t): 
K{t) = C{t)/mkBT where T is the temperature which 
characterizes the heat bath and fcs is the Boltzmann con- 
stant. The introduction of memory friction changes the 
shape of the velocity autocorrelation function consider- 
ably: it is no longer restricted to the exponentials. 



We wish to demonstrate how the process described in 
the preceded Sections can be applied as a model for the 
driving stochastic force in Ea. ljKil) . For that purpose we 
choose ACF possessing the tail ^ 1/t which characterizes 
e.g. the noise-induced Stark broadening ^23i | and nuclear 
collisions in the framework of a dynamical model |24j . 
It can be also found in problems connected with phe- 
nomena in disordered media This form of ACF is 
of special importance for molecular dynamics because it 
corresponds to the problem of scattering inside a peri- 
odic lattice ,2^]. Let us then consider the ACF given by 
Ea. (|14|l . Moreover we assume (F(t)) = 0. In this case 
K{s) = ln(l -I- 8/s)/8 and Eq.lUHl) reads 



R{s) 



1 



s + ln(l-h8/s)/8 



(19) 



In order to obtain the resolvent R{t) we need to inverse 
the above transform. Computing the usual contour inte- 
gral produces the following result 



R{t) 



— at , 
8 



e (ci sin ht + C2 cos ht)- 



dx 



(20) 



5a;-ln(8/a:-l)]2 + ' 



where the constants a = 0.3511, b = 0.2995, ci = 0.2297 
and C2 = 1.603 follow from the numerical evaluation of 
poles in the Ea. H19|) . The resolvent R has the interpre- 
tation of the velocity autocorrelation function C^,, 



C,{t) - (v(O)v(t)) = ^R{t) 



(21) 



R(t) falls from i?(0) = 1 to negative values and then 
rises, approaching zero very slowly from below. The 
behaviour of Cy(t) at t — > oo is determined by the in- 
tegral in Ea. (|20|) . In this limit, it becomes simpler: 
^ e^*^/ In^ a; dx. Integrating over t yields the inte- 
grand in the form e~*^ /{x\n' x) and the integral over x 
can be estimated |0 as ~ 1/lnt. The final expression 
reads 



-1 
tint 



{t oo). 



(22) 



Therefore the tail of Cy (t) diminishes very slowly, like the 
tail of C{t), and it is negative. 

The velocity autocorrelation function determines the 
transport properties of the system: the diffusion coeffi- 
cient can be expressed in terms of the Laplace transform 
of Cy{t) in the form: V — C.i,{s — 0). Since for C{t) given 
by Ea. ((n|l V — 0, the transport is subdiffusive. We 
come to the same conclusion by the direct calculation of 
the position variance {r'^){t). Integrating Cy{t) twice over 
time, we get the following estimation: 



(23) 



Therefore the subdiffusion is very weak and hardly dis- 
tinguishable from the normal diffusion. The same form 



5 



of the anomalous diffusion has been found in a chaotic 
(deterministic) system 26] and it has been attributed to 
the intermittency. 

Our aim is to study the motion of the particle by a di- 
rect simulation of stochastic trajectories, assuming that 
the driving force in Eq. H16|l is modeled by means of the 
stochastic process x(t) and satisfies Eq.(^. We restrict 
our analysis to the case V{r) — 0. Inserting of the so- 
lution of Eq. into Eq. (|17|l yields the two-dimensional 
trajectory of the particle velocity: 



v(i) = R{t)^r{0)+ 

r rt-t,^ " 

-m^^ x„+i Xfc 



fc=l 



R{T)dT , 

(24) 



where by sampling of the consecutive jumping times tk 
we apply Eq.Q with a = 0.5. Moreover, in the following 
we take the kernel width a — 2.5. A simple quantity one 
can evaluate from the Ea. (|24|l is the time dependence of 
the velocity variance (v^(t)) where the average is taken 
over the stationary distribution of the random force 
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FIG. 2; The velocity variance calculated from Eg. 125^ (solid 
line) and by numerical simulation from Eg. 12411 (dots). We 
assumed: v(0) = 0, T = 1 and m = 1. 

Fig. 2 presents this quantity, calculated with the initial 
condition v(0) = 0, for T = 1 and m — 1. On the other 
hand, the velocity variance can be derived analytically 
from the Ea. H24|) : the expression for (v^(t)) involves only 
the second moment of the noise: 



■ /* f R{T)R{T')C{\T-T'\)dTdT'. 

Jo Jo 



(25) 



The velocity variance appears to be independent of a spe- 
cific noise model and then analytical and numerical re- 
sults should coincide. Indeed, Fig. 2 demonstrates very 
good agreement of both results; they indicate the relax- 
ation to the thermal equilibrium (v^) = fc^T/m which is 
apparently reached at about t = 4 |23 • 
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FIG. 3: Time evolution of the probability density distribution 
of the first velocity component v^. The stochastic ensemble 
consists of 5 ■ lO" trajectories for each time. 



In a similar way, utilizing Eg. 124(1 . we can determine 
a density distribution p(v, t) which means a probabil- 
ity that the velocity of the Brownian particle is in the 
interval (v, v + dv). Fig. 3 presents this distribution, cor- 
responding to the first velocity component v^, for large 
times. The central part of the distribution is equilibrated 
already at i = 5 but tails are not yet developed; they ter- 
minate with high and narrow peaks which originate from 
the initial condition p(vx, 0) — S(vx). At short times (not 
shown in the figure) the peaks are still higher and expand 
gradually with time from a vicinity of the point Vx — 0- 
Full relaxation of the tails - which fall off faster than 
the Gaussian - to the stationary distribution is achieved 
at i = 20. Nevertheless, the memory about the initial 
condition is preserved for a very long time. The distri- 
bution of the second velocity component Vy, presented in 
Fig. 4, looks different; the width is much smaller and the 
tails show the exponential shape. A complete relaxation 
to stationary distribution is reached already at t = 10. 
The difference between the distributions for both veloc- 
ity components follows from anisotropy of the function 
1/(0): there are no infinite waiting times corresponding 
to the motion in the y direction. 

The energy spectrum of the Brownian particles devi- 
ate considerably from the Maxwellian distribution. Fig. 5 
presents the time evolution of the probability density 
distribution of the energy E — 0.5(ti^ ~^ "^y)- small 
values of the energy the curves have a cusp, whereas 
the tail of the distribution corresponding to the equi- 
librium state can be parameterized by the function 
0.5 exp(-0.5£'^) {E > 2). It is interesting that the prob- 
ability density function which characterizes the transport 
dynamics in the framework of the continuous-time ran- 
dom walk, predicts a similar cusp for the subdiffusive 
motion 0. 
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FIG. 4: The same as Fig.3 except for the second velocity 
component Vy. 
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FIG. 5: Time evolution of the probability density distribution 
of the energy E = 0.5(v^ + Vy). 

V. SUMMARY AND DISCUSSION 

The jumping process presented in this paper is char- 
acterized by the jump size probability distribution and 
the waiting time distribution which are correlated. The 



jumping rate depends on the process value which is 
kept constant between consecutive jumps. The process 
is Markovian and stationary, the corresponding master 
equation possesses a nontrivial time-independent solution 
which is completely determined by the jumping rate and 
which docs not depend on the jump size distribution. We 
have studied the process in its two-dimensional version 
for jumps which do not change the norm of the process 
value. The expression for ACF with power-law tails has 
been derived. We have demonstrated that it is possible to 
construct in a simple way a process which is a 1//— like 
noise. 

Despite the fact that the waiting time distribution is ex- 
ponential, the intervals of constant process values can 
be very long and actually algebraically distributed. This 
conclusion is not surprising because the mean value of 
the exponential distribution is also a stochastic variable. 
Then the existence of long tails of the waiting time distri- 
bution does not rule out a relaxation to the equilibrium. 
The procedure described in the paper allows us to con- 
struct stochastic trajectories corresponding to a wide 
class of power law ACF in a simple manner. Therefore 
it can serve as a model of physical phenomena and can 
be used as a stochastic force in the generalized Langevin 
equation. We have solved this equation for an exemplary 
form of ACF, ~ utilizing our process. Since waiting 
times are correlated with the direction of the noise vec- 
tor, the resulting velocity distribution exhibits a strong 
anisotropy. The distribution of the first component, cor- 
responding to long waiting times, has rapidly falling tails 
and indicates an extremely long memory about the ini- 
tial condition, despite the fact that the comprehensive 
shape of the distribution equilibrates relatively fast. On 
the other hand, the tails of the distribution correspond- 
ing to the second component coincide with the standard 
Gaussian. 

The tail of ACF is determined predominantly by the long 
waiting times and then only one component of the pro- 
cess value is crucial for its shape. Therefore, this com- 
ponent can constitute a one-dimensional counterpart of 
our two-dimensional jmnping process which still has the 
power-law ACF. This remark is important if one requires 
a model of noise possessing an arbitrary dimensionality. 
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